function [vf,vr,a]=steadystate(kf,kr,Sf,Sr,epsilon,setpoint,num_fluxes,num_metabolites)
% generic code intended to work with any network
% this code is not optimized for paralellization and will run much slower
% than needed. The same structure could be used to generate ODEs instead of
% calculating metabolites concentrations using the full stochiometric
% matrices for forward and reverse reactions
% kf and kr are forward rate constants
% Sf and Sr the substrate stoichiometric matrices
% epsilon is a small constant to prevent division by zero, and can also be 
% used to moderate the amount of scaling by "a"

%t=1; timestep not implemented yet -> t=1
n= num_fluxes;
p= num_metabolites;

a=ones(1,p);
vf=ones(1,n);
vr=ones(1,n);
m_old=ones(1,p);
m=10.*ones(1,p);

% This routine is run untill steady state is reached. This is determined by
% reaching a setpoint value of the lowest "a" value of all metabolites.
% Because "a" is inversely related to metabolite concentraion change, one
% could either use max(m) or min(a)
while min(a)<setpoint
    for i=1:n % cycle through all fluxes
        vf(i)=kf(i);
        vr(i)=kr(i);
            % take the product of all metabolite concentrations raised to
            % their stoichiometries. Because the matrix is sparse, most
            % stoichiometries are 0, resulting mainly in inconsequential 
            % multiplications of 1
            for j=1:p
                vf(i)=vf(i)*m(j)^Sf(j,i);
                vr(i)=vr(i)*m(j)^Sr(j,i);
            end
    end
    for j=1:p % cycle through all metabolites
        m_old(j)=m(j); % store old values to calculate relative change in metabolite concentrations
        for i=1:n % cycle through all fluxes 
            % add/remove all fluxes to metabolites participatig in a reaction
            % the flux is multiplied by the stoichiometry
            m(j)=m(j)-vf(i)*Sf(j,i)+vr(i)*Sr(j,i); 
        end
        
        a(j)=1/(epsilon+abs(m_old(j)-m(j)));
    end
end

    
                    
        